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Thermodynamics of a spin-1 Bose gas with fixed magnetization 
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We investigate the thermodynamics of a spin-1 Bose gas with fixed magnetization including the 
quadratic Zeeman energy shift. Our calculations are based on the grand canonical description for the 
ideal gas and the classical fields approximation for atoms with ferromagnetic and antiferromagnetic 
interactions. We confirm the occurence of a double phase transition in the system that takes place 
due to two global constraints. We show analytically for the ideal gas how critical temperatures 
and condensed fractions are changed by a non-zero magnetic field. The interaction strongly affects 
the condensate scenario below the second critical temperature. The effect imposed by interaction 
energies becomes diminished in high magnetic fields where condensation, of both ferromagnetic and 
antiferromagnetic atoms, agree with the ideal gas results. 
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I. INTRODUCTION 

A spinor Bose-Einstein condensate (BEC) is a multi- 
component condensate with an additional spin degree of 
freedom, which has provided exciting opportunities to 
study experimentally quantum magnetism, superfluidity, 
strong correlations, coherent spin-mixing dynamics, spin- 
nematic squeezing, entanglement etc, most of them in 
non-equilibrium situations (see 7] |g|]). Despite success¬ 
ful experimental developments on spinor Bose-Einstein 
condensates, our knowledge remains limited regarding 
equilibrium properties and in particular the thermody¬ 
namics of such a gas. The main reason is a long time 
needed to reach an equilibrium state, typically several 
seconds or tens of seconds, that may exceed the lifetime 
of the condensate [?j. Nevertheless, recent experimen¬ 
tal developments allowed for investigation of the ground 
state of an antiferromagnetic spinor condensate opening 
the paths to study in details its properties at thermal 
equilibrium [§]. 

The condensation of atoms with total spin F = 1 
trapped in the three hyperfine states tuf = 1,0,—1 in 
the absence of magnetic field was investigated theoreti¬ 
cally for the first time by Isoshima et al. j7f|. The dou¬ 
ble condensation phenomenon was predicted in the pres¬ 
ence of two global conserved quantities: the total num¬ 
ber of atoms N and the magnetization M. A conden¬ 
sate starts to appear in the highest mj? = 1 component 
for temperatures below the first critical temperature and 
simultaneously in the rest two components for tempera¬ 
tures below the second critical temperature. Analytical 
expressions for the two critical temperatures and con¬ 
densate fractions were given for the ideal gas and zero 
magnetic field SE3- The condensation of interacting 
spin-1 Bose gas was considered numerically within the 
Bogoliubov-Popov approximation [9] and the Hartree- 
Fock-Popov approximation [ll|. In the latter, authors 
confirmed the double phase transition for antiferromag¬ 
netic interactions, but found a more complicated phase 
diagram for ferromagnetic interactions with a possible 
triple condensation scenario. The only one experimental 


work of Pasquiou [l2j touches the problem of the thermo¬ 
dynamics in chromium atoms with total spin F = 3 but 
for free magnetization. Indeed, for low magnetic fields 
when the magnetization is approximately conserved the 
experimental results confirm the occurrence of a double 
condensation. 

In this paper, we reconsider the topic of condensation 
in the system of spin-1 bosons with fixed magnetization. 
The ultra-cold gases are almost perfectly isolated in the 
experiment and conservation of magnetization plays a 
major role. The magnetic dipole-dipole interactions, that 
may change the magnetization, are realtively weak and 
can be neglected for F = 1 sodium or rubidium spinor 
Bose-Einstein condensates. 

We examine the thermodynamics of the ideal gas in 
the presence of quadratic Zeeman effect within the grand 
canonical ensemble. A non-zero magnetic field introduces 
a new phase in the phase diagram of critical tempera¬ 
tures that we characterize by the threshold magnetiza¬ 
tion. The condensation scenario predicted by Isoshima 
is present for magnetizations larger than the thereshold 
magnetization. When the magnetization of the system 
is smaller than the threshold magnetization, atoms start 
condensing first not in the highest tuf = 1 component, as 
it was the case for zero magnetic field, but in the mj? =0 
component. That trivial effect is present due to the shift 
of the lowest energy level of the mp = 0 component be¬ 
low the lowest energy level of the rriF = 1 component. 
We give an explicit expression for the threshold magne¬ 
tization. 

We study the interacting gas within the classical fields 
approximation [l2j combined with the Metropolis algo¬ 
rithm Q The method was successfully used to in¬ 
vestigate thermal effects in the single-component Bose- 
Einstein condensates including thermodynamics m 
vortex-dynamics (iff, critical temperature shift 17j, 
spin-squeezing fl8l |. solitons or Kibble-Zurek mecha¬ 
nism [19j, and many others, some of them reviewed 
in [20j. This numerical method includes all non-linear 
terms present in the Hamiltonian at the expense of in¬ 
troducing a free parameter that has to be well chosen. 
In this paper we explain how to adapt the Metropolis 
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algorithm for a spin-1 gas with fixed magnetization. To 
demonstrate the validity of the proposed algorithm we 
compared results of simulations with exact results for 
the ideal gas and with the approximated Bogoliubov the¬ 
ory for antiferomagnetic interactions. We confirmed dou¬ 
ble condensation for both ferromagnetic and antiferro¬ 
magnetic interactions. The condensation strongly differs 
from the result for ideal gas below the second critical 
temperature. In the high magnetic field limit, when the 
quadratic Zeeman energy dominates over the interaction 
energy, details of condensation do not depend on the in¬ 
teraction sign and are well described by the ideal gas 
results. 


II. THE MODEL 

We consider a dilute and homogeneous spin-1 Bose gas 
in a magnetic field. We start with the Hamiltonian H = 
Hq + Ha, where the symmetric (spin-independent) part 
is 


H 0 = 



i’j- (!) 


Here the subscripts j = — ,0, + denote sublevels with 
magnetic quantum numbers along the magnetic field axis 
trip = —1,0,+1, to is the atomic mass, n = ^rij = 
i’jipj is the total atom density. The spin-dependent 
part can be written as 


H a 


fd 3 r 

J2 E jnj + |:F 2 : 
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( 2 ) 


For a sufficiently weak magnetic field we can approxi¬ 
mate Zeeman energy levels by a positive energy shift of 
the tof = ±1 sublevels 5 = (E+ + E- — 2Eq)/2 ss qh 2 , 
where h is the magnetic field strength and q = (<?/ + 
9 J) FbI^E hfs> gj an d gi are the gyromagnetic ratios 
of the electron and the nucleus, /.ib is the Bohr magneton, 
£hfs is the hyperfine energy splitting at zero magnetic 
field. Finally, the spin-dependent Hamiltonian @ be¬ 
comes 


Ha 



qh 2 (n+ + n-) + y : F 2 : 


(5) 


where F 2 = n+ — n_ and F]_ = 2\ip + ipQ + ip o'*/’— | 2 are 
the square of magnetization density and the square of 
transverse spin density, respectively. In spinor conden¬ 
sates realized in laboratories, ao and 02 scattering lengths 
have similar magnitude. The spin-dependent interaction 
coefficient C 2 is, therefore, much smaller than its spin- 
independent counterpart cq. For 23 Na condensate their 
ratio is about 1:30 and is positive (antiferromagnetic or¬ 
der), while for 87 Rb condensate it is 1:220 and is negative 
(ferromagnetic order). 

By comparing the kinetic energy with the interac¬ 
tion energy, we can define the healing length £ = 
2nh/y/2mcon and the spin healing length = 

2nh/y/2mc2n. These quantities give the length scales of 
spatial variations in the condensate profile induced by the 
spin-independent or spin-dependent interactions. Here 
we consider system sizes smaller than the spin healing 
length in order to avoid a domain formation. A good 
basis for such a homogeneous system is the plane wave 
basis. 


III. THE IDEAL GAS 


where Ej are Zeeman energy levels, F = 
W’Vad/W M,^fz1p) T is the spin density, fx,y,z 
are spin-1 matrices, ip = (ip+, ipo, ip-) T : and :: de¬ 
notes the normal order. The spin-independent and 
spin-dependent interaction coefficients are given by 
Co = 4:Trh 2 (ao + 202 )/3m and C 2 = ATth 2 (a ,2 — do)/3m 
respectively, where as is the s-wave scattering length for 
colliding atoms with total spin S. The total number of 
atoms 



and the magnetization 

M = J ( n + — n-)d 3 r (4) 

are conserved quantities. 

The linear part of the Zeeman shifts Ej induces a ho¬ 
mogeneous rotation of the spin vector around the direc¬ 
tion of the magnetic field. Since the Hamiltonian is in¬ 
variant with respect to such spin rotations, we consider 
only the effect of the quadratic Zeeman shift. 


We consider a uniform gas of non-interacting atoms 
(co = C 2 = 0) with hyperfine spin F = 1 in a homo¬ 
geneous magnetic field h within the grand canonical en¬ 
semble, taking into account the quadratic Zeeman effect. 
The effective Hamiltonian of the system is 

H ef f = ^2 ^2{ e k + m 2 F qh 2 ) n k , mF ~ h N ~ , 

mp-1,0,-1 k 

( 6 ) 

with 

N = N + + N 0 + N_ 

mF k 

M = N + - 7V_ = Y2 m F n ^,rn F ■ (8) 

mp k 

Here k = 2tt/L( n Xl n v ,n z ), L is the system size and 
ni = 0,±1,±2... are integers, n k , m F are occupation 
numbers of atoms of energy e k = h 2 k 2 /2m. The chemical 
potential /i and the linear Zeeman shift ij are Lagrange 
multipliers enforcing the desired total atom number N 
and the magnetization M respectively. N mp is the num¬ 
ber of atoms in the tof th component. We consider a 
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positive magnetization M > 0 and a positive Zeeman 
energy shift qh 2 > 0 . 

The non-zero magnetic field removes the degeneracy of 
energy spectra: 


A. Transition temperatures and condensate 
fractions 

1. For qh 2 < r/ 


£k,+ = ek - n~r) + qh 2 , 
£k,o = £k — /x, 

-Ek - = £k-/^ + r? + gE 2 . 


(9) The first phase transition occurs for z+ —> 1 (or p, —^ 
( 10 ) qh 2 — 77 ) when the tuf = 1 component starts condensing. 
That Nfi 3> 1 can be seen from (1181) . The number of 
thermal atoms is then 


The ratio between 77 and qh 2 determines the order of 
energy levels. The lowest energy level is E+ for qh 2 < rj, 
or Eq for qh 2 > 77 . In addition, two effects determine the 
state of the system: (i) the occupation number imbalance 
enforced by the fixed magnetization N + > iV_, and (ii) 
the ground state energy level (Eq or E+) controlled by 
the magnetic field. 

Since the Hamiltonian is diagonal, we may calculate 
the grand canonical partition function 

S = e -P E k,m F nk,m F ( ( 12 ) 

m F ,n ktmF 

where /3 = l/fc^T, T is the temperature and Eg is the 
Boltzmann constant. The ensemble average of the occu¬ 
pation number n^ mF is 


^k,m p — 


1 <91n£ 

0 dE KmF ’ 


which gives 


^k,mp 


y p /^ e k 

1 - Zmpe-P^ 


with effective fugacities 


, _ e P^+v-qh 2 ), 

zo = , 

z _ = eP^-v-qb 2 ) _ 


(13) 


(14) 


(15) 

(16) 
(17) 


In the thermodynamic limit, keeping only dominant 
terms O(N), expressed in terms of fugacities, the number 
of atoms in the lowest energy level of each tuf component 
is 


N rn F = T~TV — ’ ( 18 ) 

-L Zm p 

while the number of thermal atoms in each tuf compo¬ 
nent is 


Nmp ~ 9§( z m F ) j (19) 

where A dB = h/\Z2irmkBT is the de Broglie wave length, 
and (jj(x) = XivS xU jn? is the Bose function. 



( 20 ) 

( 21 ) 

( 22 ) 


with z v = e |8r h The first critical temperature T c 1 can 
be obtained from the following equations: 



where we have introduced the notation z vc i = z v (T c \) 
and 


Ft(T, z v ) = <73 (1) + 33 (e^\) + <7| (z 2 ) . (25) 

Below T c 1 , only the mp = 1 component condenses. It 
is justified to assume N+ ~ N c . Then relation N c = N — 
Nj defines the condensate fraction of the tuf = 1 
component 


jV£ / T \5 Ft{T,z v ) 

~n~ 1 ~\tTi) f+ 2 (t c1 ,z vc1 ) ■ 
2 


(26) 


The second phase transition occurs for z v —> e~^ qh 
(77 —v qh 2 ) when Nq 1 and JVf. —> e~ 2 P qh /(I — 
e -2/3qh ^ j n regime, T < T c 2 , thermal populations 
are 

(27) 

jV » t = (^)’ 9 « (1) ' (28> 

(29) 

The second transition temperature T c 2 can be obtained 
using the difference between the total atom number 
N and the magnetization M. For temperatures T £ 
[T c 2 ,T c i], in the absence of condensates in tuf — 0,-1 
components, the difference is N — M ~ 2 Nfi + Nq. The 









4 


second transition temperature expressed in terms of Bose 
functions present in equations m and (1^1) is 


ksTc 2 


2irh 2 N -M 

mL 2 Ga(T C 2 ) 
2 K 7 


(30) 


where 

G i (T) = g i (l) + 2 9 3 (e- 2 ^ 2 ) . (31) 


Below T c 2 , the Bose-Einstein condensate can be formed 
in all components and condensate fractions satisfy the set 
of equations: 


N c = 



M C {T) + {N-M) 
N± = M C (T ), 



Ga(T)' 

GaVcs) 

7 J 


N't 


+ 


-2pqh 2 

w 


-2sinh(/3g/i 2 )e- /3 « /l2 . 


(32) 

Here we have introduced the condensate part of the mag¬ 
netization M c = M — M t and the thermal part of mag¬ 
netization 


M t (T) 



5| (!)-<?! (e 2 ^ 2 ) 


(33) 


A derivation of eqs. (l32l) is included in Appendix [A] An 
analytical solution of eqs. (15^1) is presented in Appendix 
m We have checked the validity of the analytical solution 
against the self-consistent numerical result. 


2. For qh 2 > g 

First, one should notice that this case does not exist 
in the absence of an external magnetic field, since rj can 
take positive values for M > 0. That is a new area of 
the phase diagram, which appears due to the quadratic 
Zeeman effect. 

The first phase transition. This time, the tuf = 0 
component undergoes condensation first, which means 
that zq —!> 1 (or p —> 0) and N§ 3> 1. One obtains new 
expressions for NT, Nq and /W, which hold under the 
critical temperature T c i: 

Nl = (^) g^e-^ 2 z v ~ l ), (34) 

N ° T= ( a ^) (35) 

^=(^) 3 <?f ( e-^\). (36) 

The first critical temperature T cl and the fugacity at the 
critical point z r)c \ can be obtained from the following 


equations: 


( h 

\^dB{T c l) 

( L 

\^dB(T c l) 


Fl{T cU z vcX ) , (37) 

ff|(^ci _ 1 e T fc')-ga^cie 1 ^) , 

(38) 


where we have introduced 

F|(T, z v ) = g^z^e-^) + gz(l) + g^e-^ 2 ). 

(39) 

At the critical the point the fugacity is smaller than one 
{z v (T c i) < 1) since M > 0 and g3 is an increasing func¬ 
tion of its argument and takes positive values. 

Assuming that Nq ~ N c for T £ [T c2 ,T c i], once again 
the relation N c = N — VA Nff defines the condensate 
fraction in the tuf = 0 component 


Nq ( T \* F U T ’ Z ri) 

F 2 (T cl ,Zr, cl ) ' (40) 

2 

The second phase transition. One expects z v ~ . 

implying that the tuf = 1 component starts condensing: 
N% > 1 and again 7V£ e~ 2/3qh2 /{I - e~ 2f3qh2 ). Never¬ 

theless, in this regime we should define T c 2 in the other 
way. Neither N nor N — M can be used anymore, since 
they involve Ay /N which is now unknown in the interme¬ 
diate region of temperatures. The only solution is to use 
the magnetization M , and define T c 2 as the temperature 
for which Nf ~ N± « N, that is 


ksT c 2 — 


27 xh 2 
mL 2 



M 


— 2 qh,2 

g3 (e 



(41) 


what is equivalent to M = M t {T C 2 ). This choice 
is justified in the thermodynamic limit when N » 
e - 2 pqh — e -2/3g/i ^ anc j mathematically within our 
equations for f3qh 2 1 and any N. 

Below T c2 , condensate fractions satisfy the set of equa¬ 
tions: 


N c = N - 
jV|- ; 
2 _ j_ 
NS NS - 


(lV 
\^dB J 
= M C (T) 

-2Pqh 2 


N° 


2fi , f(l) + g\ ( e 2/3qh2 ) 

5 

-2sinh(/J g /i 2 )e- /3 ® ,t2 


(42) 


Different from (15121) is the first equation of (1421) only. 


B. Phase diagram 

The non-zero magnetic field changes dramatically the 
phase diagram of the critical temperatures, which is 
shown in Fig. [0 The phase diagram consists of four 
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FIG. 1. (Color online) (a) and (b) show the phase diagram of the critical temperatures. T c i is marked by solid lines and 
T c 2 is marked by dashed lines. In (a) black lines are for qh 2 = 0, red lines are for qh 2 = 0.5h 2 /mL 2 and blue lines are for 
qh 2 = 12.5 h 2 /mL 2 . T® = (2nh 2 /mL 2 )(N/((3/2)) 2 ^ 3 is the critical temperature for one component condensate in the box 
potential. In (b) the same for qh 2 = 12.5fi 2 /2mL 2 . The arrow indicates the threshold at the critical temperatures intersection 
point. Particular parts of the diagram are "/I"-thermal atoms (no condensate), "^''-condensate in the component tof = 1, 
"S'"-condensate in the component tof = 0 only, "C"-condensate possible in all components, (c) The threshold temperature 
Tth/Tc (solid line) and the threshold magnetization M t h/N (dash-dotted line) at the critical temperatures intersection point 
as a function of the quadratic Zeemann energy shift qh 2 /e with e = fi 2 /2mi 2 . The asymptotic values of the threshold critical 
temperature T////Tf h = 2~ 2 + and the threshold magnetization M/f/ = 1/2 are marked by dashed and dotted lines respectively. 
Here, the total number of atoms is TV = 10 4 . 


phases ” A” , ” B v , ” B'” , ” C” separated by the two criti¬ 
cal temperatures T c \ and T c 2 . Depending on the value 
of the temperature, the system can be: ”A”- a non¬ 
degenerate thermal gas, ” B”- condensate in the compo¬ 
nent mj? = 1 or ’’.^’’-condensate in the rriF = 0 and 
thermal atoms in other components, ”C”- a condensate 
in m f = 0 and mp = 1 while in the mf = — 1 is a gas 
with non-negligible fraction of atoms in the lowest energy 
level for fiqh 2 C 1, or with negligible fraction of atoms 
in the lowest energy level for f3qh 2 1. A possible desti¬ 
nation is controlled by the magnetization, with a special 
role of the threshold magnetization at the critical temper¬ 
atures intersection point M t h = M(T = T c \ = T c 2 ). If 
M < M t h then to obtain particular quantities, one should 
use expressions from subsection IIII A 21 in the opposite 
case (M > Mth) from subsection IIII A 1~1 The procedure 
to obtain numerical values for the critical temperatures 
is explained in Appendix [Cj 

Analytical expressions for the threshold critical tem¬ 
perature T t h and the threshold magnetization M t h are 



M th 

N 


N 

2ff3/2(l) +93/2 (e~ 2qh2 ! kBTth ) 


(43) 


3ff3/ 2 (l) (Tth\ 3/2 

n y c J 


(44) 


where C = h 2 /2'KmL 2 kB- The above expressions are ob¬ 
tained by comparing critical temperatures T c \ and T c 2 
for both qh 2 > 77 and qh 2 < 77 . In fig. QJ we show the 
threshold critical temperature (1431) , the threshold magne¬ 
tization (PHI) and their asymptotic values for f3qh 2 —> 00 , 
they are T%/T° h 2 " 2 / 3 and M°£/N 1/2 respec¬ 

tively. 


C. Condensed fractions 

Condensate fractions, solutions of equations (1521) for 
M > M t h and solutions of equations PI21) for M < M t h, 
are presented in fig. [2] and [3] respectively and they are 
marked by lines. Points are results of the Metropolis 
algorithm adapted to the model (more details concerning 
the algorithm can be found in section HV Bl) . 

Figure [2] is for values of the magnetic field in the area 
M > M t h where qh 2 < 77 . These graphs show that 
modifications of condensed fractions occur mainly for 
low magnetic fields. The effect of the non-zero magnetic 
field is the most visible on the condensate fraction in the 
rriF = — 1 component. Notice, at zero magnetic field the 
fraction of the condensate in the mp = — 1 component 
decreases simply with the temperature, see fig.[ 2 } 7 . In the 
transient magnetic field regime, the condensed fraction in 
the rriF = — 1 component increases from zero, reaches a 
maximum and then decreases to zero at the second crit¬ 
ical temperature, see figHJ). The condensed fraction in 
the rriF = — 1 component decreases quickly and can be 
neglected in the high magnetic field regime, see fig. ( 2 J:. 
Condensed fractions are linked together, when the con¬ 
densed fraction of the = -1 component disappears, 
in the meantime the condensed fraction in the rriF = 1 
component decreases and the condensate fraction in the 
wif = 0 increases. Nevertheless, the slope-breaking that 
occurs at T c 2 , already present when h = 0, is still neat. 
The fugacity varies dramatically near the zero temper¬ 
ature for small values of magnetic fields, what explains 
sharp variations of the condensed fractions in fig. [ 2 ]:). 

Figure [3] is for the magnetization M < Mth where 
qh 2 > 77 . The value of magnetization is M = 50 and 
is very small as compared to N = 10 4 , therefore the dif¬ 
ference between N1 and N l _ i is not visible. Notice, strong 
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FIG. 2. (Color online) Condensate fractions for M > Mth, with N = 10 4 and M = 6 x 10 3 ; (a) qh 2 = 0, (b) qh 2 = 0.1 h 2 /mL 2 
and (c) qh 2 = h 2 /mL 2 . n c = N c /N is the total condensate fraction (solid black line), is the condensate fraction in rriF = 1 
component (dashed red line), and rig, n c _ in ttif = 0 (dot-dashed green line), rriF = — 1 (dotted blue line) respectively. Lines 
are the solution of equations (I32|) while points are results of Monte Carlo simulations. 



FIG. 3. (Color online) Condensate fractions for M < Mth, with N = 10 4 and M = 50; (a) qh 2 = 0 and (b) qh 2 = h 2 /mL 2 . 
n c = N c /N is the total condensate fraction (solid black line), n+ is the condensate fraction in the ttif = 1 component (dashed 
red line), and rig, ni in rriF = 0 (dot-dashed green line), rriF = — 1 (dotted blue line) respectively. Lines are solution of 
equations (g2D while points are results of Monte Carlo simulations. Particular points in (a) correspond to averaging over 
different representations of an ensemble and show strong fluctuations of condensate fractions in the regime of zero magnetic 
field and almost zero magnetization. 


fluctuations of the condensate fractions for zero magnetic 
field that are results of Monte Carlo simulations, see dif¬ 
ferent points in fig. [3jr. Indeed, for zero magnetization 
and zero magnetic field the ground state of the ideal gas 
is strongly degenerate El , what gives rise to strong fluc¬ 
tuations of condensate fractions. The non-zero magnetic 
field reduces degeneracy and hence reduces fluctuations 
of condensate fractions in fig. [3 Jd. 

IV. THE INTERACTING GAS 

The ground state of a spin-1 Bose gas in the presence 
of ferromagnetic and antiferromagnetic interactions was 
widely studied within the single-mode aproximation [2l| 
and beyond (221, and was investigated in experiments 
for antiferromagnetic condensates J8J. The structure of 
the ground state is quite complex and depends not only 
on the magnetization and magnetic field but also on the 
relative phase between components of the Bose gas. It 
consists of a polar, nematic or magnetic state, two com¬ 
ponent or three component solutions with phase and anti¬ 


phase matching for ferromagnetic and antiferromagnetic 
interactions respectively. We are aware of the temper¬ 
ature dependence of such structures, in particular the 
boundaries between different phases. 

The non-zero temperature introduces a multi-mode 
structure, therefore we describe the system within 
the classical fields approximation that takes into ac¬ 
count thermal populations and interactions among many 
inodes. Indeed, classical fields and stochastic methods 
[23j as well as Hartree-Fock or Hartree-Fock-Popov ap¬ 
proximations [24| were applied for spinor condensates at 
non-zero temperature but for free magnetization. Among 
all of finite temperature methods that are used for single¬ 
component condensates, the classical fields-like are not 
perturbative and thus contain all nonlinear terms that 
are present in the Hamiltonian. It makes them very suit¬ 
able for study thermodynamics in the whole temperature 
range, what is not a case for methods based on the Bo- 
goliubov approximation. 

Below we just briefly remind the main concept of the 
classical field approach, more details concerning Hie foun¬ 
dations of the approximation can be found in [l3|. 
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A. The classical fields approximation 

The classical fields approach consists of (i) replacement 
of the creation and annihilation operators by complex 
amplitudes, (ii) restriction of the summation over modes 
to a finite number extended all the way to the momentum 
cut-off K max . 

The field operator is replaced by a classical field (com¬ 
plex function) of well-defined number of momenta modes: 

</b( r ) = E a i( k )- ( 45 ) 

k<K ma * 

The energy E^ of such a classical field is given by dis¬ 
cretization of the hamiltonian H = Hq + Ha, eqs. (CD 
and ©. The total number of atoms is 

N = J2 E M k )| 2 (46) 

j k<K max 

and the magnetization is 

M= J2 (Mk)| 2 -|a-(k)| 2 ). (47) 

k<K max 

Various observables have a more or less pronounced 
dependence on the cut-off. Here we choose the cut¬ 
off momentum such, that in the thermodynamic limit 
the non-condensed density for a single component ideal 
Bose gas in degenerate regime is exactly reproduced 
by the classical field model [25|. The condition gives 
EjA max ~ 2.695 ksT, where EjA max = 3h 2 (ir / L) 2 /2m is 
the maximal kinetic energy on the grid. 


B. The Metropolis algorithm for a spin-1 Bose gas 
with fixed magnetization 

We adapt the Metropolis sche me [bill to the system 
of classical fields as described in [15j]. The main idea 
of this Monte Carlo method is to generate a Markovian 
process of a random walk in phase space. All states of the 
system visited during this walk become members of the 
statistical ensemble and are used in ensemble averages. 

In order to obtain a statistical average of any observ¬ 
able A: 

A i = ^Y.{4 ] \M s) ) ( 48 ) 

s— 1 

one should generate J\f copies of the classical fields . 
A canonical average is obtained in the limit of Af —> 
oo provided the number of members of the ensemble 
with energy E ^ is proportional to the Boltzmann fac¬ 
tor e ~ E */ k BT _ This can be achieved in a random walk 
where a single step of the Markov process is defined as 
follows: 


1. A set of amplitudes (k) determines the state se¬ 
lected to be a member of the canonical ensemble at 
the sth step of the random walk. The correspond¬ 
ing energy E^ of the classical field is calculated 
according to © and ©• As the initial condition 
(s = 1), any state that satisfies the condition of 
the fixed total number of atoms N and the mag¬ 
netization M may be chosen as a member of the 
ensemble. 

2. A trial set of amplitudes a^(k) is generated by a 

random disturbance of a^(k) = a^(k) + dj s ^(k) 

followed by normalization to account the total num- 

~/ s \ 

ber of atoms. This way a trial classical field if)) ' 
is obtained. The corresponding magnetization M s , 
energy E^, the energy difference A s = E^ - E^, 
as well as the Boltzmann factor p s = e - A Ak B T are 
then calculated. 

3. If the magnetization M s of a trial set of amplitudes 
satisfy |M — M s \ < SM then a^(k) can be consid¬ 
ered as a new member of the ensemble. 

4. A new member of the Markov chain aj s+1) (k) is se¬ 
lected according to the following prescription: ( i ) if 
A s < 0 then the trial state becomes a new member 
of the ensemble a^ s+1 '*(k) = a^^k), (ii) if A s > 0 
then a random number 0 < u < 1 is generated. If 
u < p s then the trial state becomes a new member 
of the ensemble. In the opposite case u > p s , the 
"initial" state a[^ (k) is once more included in the 
ensemble a^ s+1 '(k) = a^(k) . 

The convergence of the procedure is the fastest when ap¬ 
proximately every second trial state becomes a member of 
the ensemble. This factor depends on the assumed max¬ 
imal value of displacements S^' 1 (k) which can be mod¬ 
ified during the walk. The parameter SM should be 
small enough to ensure almost constant magnetization 
M. Note, some number of initial members of the ensem¬ 
ble should be ignored in order to avoid an influence of 
the arbitrarily selected initial state of the system. 

In order to demonstrate the validity of the algorithm 
we compare Monte Carlo simulations with the exact so¬ 
lutions for the ideal gas in figures [5] and 0 and with 
the approximated Bogoliubov theory for antiferromag¬ 
netic condensate in figlU In the latter case, analytical 
solutions are given by the Bogoliubov transformation for 
antiferromagnetic interactions and are valid in the low 
temperature limit below the critical magnetic field [26|. 
Both comparisons are satisfactory what allows to use the 
algorithm in the wider range of interactions. 

C. Numerical results 

In figures [5] and [U] we show results of numerical sim¬ 
ulations using the Metropolis algorithm. Figure [5] is for 




FIG. 4. (Color online) A test of the Metropolis algorithm for 
23 Na spinor condensate in the low temperature limit. The 
condensed fractions rij, for j = dh, 0, are plotted in the figure 
as a function of relative magnetization M/N for (a) qh 2 = 0 
and (b) qh 2 = 0.01ft 2 /mL 2 . Particular colors denote conden¬ 
sate fractions in tuf = 1 (red), mp = 0 (green) and mp = — 1 
(blue) component. Solid lines: the Bogoliubov theory, Points: 
Monte Carlo results. The total number of atoms is N = 10 s . 


the magnetization M = N/2, and figure [6] for M = 50. 
Condensate fractions for atoms with antiferromagnetic 
interactions are marked by filled points while for atoms 
with ferromagnetic interactions by open ones. Particu¬ 
lar colored symbols denote condensate fractions in the 
nif = 1 component (red circles), rriF = 0 component 
(green squares) and in the rriF = — 1 component (blue 
diamonds). Solid lines denote results for the ideal gas, 
added in the figures for comparison. 

It clearly reveals the double phase transition that oc¬ 
curs in the system as it is determined by the ideal gas 
calculations. It does not seem that critical tempera¬ 
tures were affected very much by interactions. More¬ 
over, even condensate fractions for the range of temper¬ 
atures T £ [T c 2 ,T c1 ] and any magnetic field follow the 
ideal gas prediction. It is not very suprising since the 
system condenses in this regime like the single compo¬ 
nent gas. Below the second critical temperature the con¬ 
densate scenario results from the competition between 
spin-dependent interactions (dominant at low magnetic 
fields) and the quadratic Zeeman energy (dominant at 
large magnetic fields). The impact imposed by interac¬ 
tions is the most visible in the low magnetic field regime 
where ferromagnetic atoms condense differently than an- 
tiferromagnetic, and both do not match the ideal gas 
curve, see figs. [5jr and [Ur. Nevertheless, dissimilarity in 
populations of a given component between both interac¬ 
tion types is not so large. 

The antiferromagnetic interaction reduces the conden¬ 
sate population in the rriF = 0 component in all tem¬ 
perature range for magnetic fields below its critical value 
known from the ground state analysis [2lJ, see fig. [5ji. 
Simultaneously, the condensate fraction in the rriF = ±1 
components decreases like (T/T 3 ) 3 / 2 . The ferromagnetic 
interaction allows for condensation in all components, 
and populations in the lowest momentum mode may de¬ 
crease or increase up to the second critical temperature 
depending onraj?. In the other parameters regime con¬ 
densate fractions may not simply decay with the tem¬ 


perature but may also increase up to some temperature, 
reach a maximum and then decrease, see filled red points 
in fig. [5)3 for example. This feature is also observed for 
the ideal gas. In the high magnetic field regime, where 
the quadratic Zeeman energy dominates over the spin- 
dependent interaction energy, the condensate scenario 
matches the ideal gas prediction for both types of in¬ 
teractions, what can be seen in figures [5jc and (5}). 

The interesting case of almost zero magnetization and 
zero magnetic field is presented in fig. [HJi. We observe 
strong fluctuations of condensed fractions for atoms with 
antiferromagnetic interactions (shown in the inset), what 
is not the case for ferromagnetic atoms (shown in the 
main window). In the inset of fig. (Ur numerous points 
are obtained by averaging over different representations 
of ensemble members. Results of the Monte Carlo simu¬ 
lations strongly fluctuate, and additionally they are sen¬ 
sitive to the parameters of simulations (members of en¬ 
semble, or SM for example). Similarly as for the ideal 
gas, the ground state of the antiferromagnetic condensate 
is degenerated what gives rise to observed fluctuations. 
The phenomenon that is behind this effect is called spin 
fragmentation and was already investigated theoretically 
for the antiferromagnetic spinor condensate [27j]. 


V. SUMMARY 


In summary, we have studied the thermodynamics of 
a spin-1 Bose gas with fixed magnetization in the pres¬ 
ence of a non-zero magnetic field. We have given explicit 
expressions for the two critical temperatures and all con¬ 
densate fractions for the ideal gas. We have shown the 
occurence of a new phase in the phase diagram of critical 
temperatures. The interacting gas was studied within the 
classical fields approach, that is not perturbative and in¬ 
cludes all nonlinear terms present in the hamiltonian. An 
alternative method, namely stochastic Projected Gross- 
Pitaevskii equation, was lately adapted to the case of 
a spin-1 Bose gas but for free magnetization [23f. We 
find that interactions strongly affect the condensation 
scenario below the second critical temperature and for 
low magnetic fields. In this regime of parameters the 
thermodynamics of ferromagnetic and antiferromagnetic 
gases are different. The condensation is not affected 
much by interactions for values of temperatures between 
the two critical temperatures T £ [T c i,T C 2 ] for all val¬ 
ues of magnetic fields. Furthermore, the condensation is 
not affected by interactions in the whole temperatures 
range in the high magnetic field limit. Generalization to 
a Bose gas with arbitrary spin F is straightforward. Our 
results open the path to study the influence of a multi- 
mode structure on the properties of spinor condensates, 
providing an interesting direction for a future work. 
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FIG. 5. (Color online) Condensate fractions for 23 Na atoms with antiferromagnetic interactions (filled points) and for s 'Rb 
atoms with ferromagnetic interactions (open points). Particular colored symbols denote condensate fractions in the tuf = 1 
component (red circles), rriF = 0 component (green squares) and in the rriF = — 1 component (blue diamonds). Solid lines are 
results for the ideal gas. The total number of atoms is TV = 10 4 , the magnetization M = N/2 and values of magnetic fields are 
(a) qh 2 = 0, (b) qh 2 = 0.124 ti 2 /mL 2 and (c) qh 2 = h 2 /mL 2 . 



FIG. 6. (Color online) The same as in fig. [5] but for magnetization M = 50. Condensate fractions for 23 Na atoms with 
antiferromagnetic interactions are marked by filled points and for Sl Rb atoms with ferromagnetic interactions by open points. 
Particular colored symbols denote condensate fractions in the rriF = 1 component (red circles), rriF = 0 component (green 
squares) and in the rriF = — 1 component (blue diamonds). Solid lines are results for the ideal gas. The total number of atoms 
is N = 10 4 and the values of magnetic fields are (a) qh 2 = 0 and (b) qh 2 = h 2 /mL 2 . In (a) condensed fractions for atoms with 
ferromagnetic interactions are presented in the main window. Noumerous points in the inset show condensate fractions for 
atoms with antiferromagnetic interactions that are obtained by averaging over different representations of ensemble members. 
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Appendix A: Equations for condensate fractions 
when T £ [0, T c 2 ] 

Here we show how to obtain equations (1521) for con¬ 
densed fractions. 

The first one is obtained by writing N c = N1 + N£ + 
AT = N — N t in the following 



then after introducing Gj and M T it has a form 

N C = N- M t (T) - ( T ) ■ ( A2 ) 

Knowing that N — M = (L/Ads) 3 G| (T), after some al¬ 
gebra one finds the first equation of (1521) . 

The second equation of (1521) is just rewriting the total 
magnetization in terms of its condensate and thermal 
parts M = (APjl — Nt) + {N+ — Nt), and an observation 
that the whole thermal part simply reduces to M t (T). 

The third formula of (1321) is a bit more tedious to ob¬ 
tain. It is a peculiar case of a more general formula, valid 
in any regime, that we shall prove now. Starting from the 
set of equations 


(Al) 


(A3) 
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m = 


Zo 


z+z^e 


flqh 2 


0 1 — zo 1 — z+z v eP qh2 


(A4) 


N1 = 


z+z v 


1 — Z- 1 — z+z v 2 ’ 


(A5) 


one rewrites 


1 


N‘ 


+ z + 


(A6) 


N§ z+z^e^b 2 


- 1. 


(A7) 


1 


1 


N1 z+z v 2 


- 1. 


(A8) 


which shows that 


z p-Pqh -l p-pqh 


Arc 


— __i_9_ 7 P -pQ h _ z -i p-Pqh 

N c + - NS + " " 

(A9) 

and it leads to the third equation of (15^1) in the limit 
-> p-P qh2 


Appendix B: The analytical solution of equations for 
condensate fractions when T E [0, T c 2 ] 


The algebraic considerations of (15^1) lead to the follow¬ 
ing equation for Nq : 


Nf + aNf + WV n c + c = 0, 


(Bl) 


where 


1 — 2uN c + 2 u 


(B2) 


2 N c u — 2 M e ffU + 2 N c — uN c 2 + uM 2 f f 

b= ---—-----(B3) 


select the one that has the proper limit when qh 2 —> 0. 
To do it one defines: 

X = NS+1, (B5) 

which allows to put the polynomial into the form: 

X 3 +pX + q = 0 (B6) 

with 

P = b- y, (B7) 


q = ±(2a 2 -%) + c. (B8) 

Then, one writes: 

X = u + v , (B9) 

and notices that u 3 and v 3 are solutions of 


X 2 + qX - f/27 = 0 . 

Then one introduces 

A= 27g 2 +4p 3 
27 


(BIO) 


(Bll) 


Numerically, it appears that A < 0 and p < 0, which 
means that there are three solutions: 


X k = 2\J— p/3cos ( - arccos 

(B12) 

where k £ {0,1,2}. To find A}}., one has to keep in mind 
that 

N£ k =X k -^. (B13) 

Eventually, one finds that Nq x should always be selected 
for Nq because it is the only solution that gives the ap¬ 
propriate limit when qh 2 —> 0. Then, having Nq, one 
can easily calculate N+, N^_ with the first and second 
equations of (15^1) . 



Appendix C: How to obtain the phase diagram 


c = 


»C 2 - Ml,, 


(B4) 


and u(q,h,T) = sinh ((3qh 2 )e P qh - 2 , u(q,h,T) = 
cosh (f3qh 2 )e~^ qh2 . 

The quadratic Zeeman effect transformed the equation 
for Nq which was a second degree polynomial for the 
zero magnetic field, into a third degree polynomial. That 
polynomial has three roots. One needs to select, among 
those solutions, the only one that is physical: real, non¬ 
negative, and with values between 0 and N. To avoid 
numerical difficulties, one can find analytical solutions of 
this equation using for instance Cardan’s method, and 


In this appendix we explain how to compute numer¬ 
ically the transition temperatures T c i and T c 2 for any 
fixed magnetization. 

There is an additional difficulty to compute T c i com¬ 
pared to the case when h = 0, since its definition involves 
Zrjd = z v (T c 1 ), which is unknown. One can determine 
z vc 1 from the constant of motion M/N, but to do so, one 
has to know the value of T c i, as can be seen from the 
following sets of equations : 


T cl = C 


2 



(Cl) 

















11 


M 

TV 


_-ga(^ci)_ 

flf (1) + 3| (e Jcl ^ci) + g§ O^i) 


if q/i 2 3 4 5 6 7 8 9 10 < 77 , and 


Tel =C 


2 



(C2) 


(C3) 


qh 2 qh 2 

M (e fc B T <= 1 Zr,ci~ l ) - ga(e ‘b t ' 12, c1 ) 

/V qfc 2 qh 2 5 

5 |(e fc B T ‘=i^ci _1 ) +fl|(l) +3|(e fc B'3ci ^ cl ) 

( C4 ) 

if gft 2 > 77 . There are no further independent equations 
available for those two quantities, so they have to be 
solved in a self-consistent way. 

Let us suppose that the magnetization is such that the 
system is in the area where qh 2 < 77 . One knows the 
value of T c i without magnetic field, and one can sensibly 
expect that if a magnetic held is switched on, the critical 
temperature will be of the same order of magnitude as it 
used to be, so one puts the value of T c i in the absence of 
any held qh 2 = 0 to compute the value of z v c i, and then 
puts this value into (1C II) to compute the corrected value 
of T c i, that can be used in (IC2I) to compute z vc i. Those 
operations should be performed as many times as needed 
to make the effect of the wrong initial value disappear. 
The convergence is fast, and after few steps one is close 
to the hxed point for T c i. 


One can proceed in the same way using (IC3I) and (IC4I) 
if qh 2 > rj, but how is it possible to know at once if one 
is in this case or in the other? One does not know it, 
but it is of no importance whatsoever if one uses a little 
trick familiar to chemists. At the beginning one makes 
some assumption, and then checks if the computed value 
Zrjd is consistent with this guess. If not, then it means 
that the assumption was wrong, and that the system is 
in the other area. One should begin calculations again. 
Numerical problems can occur if one is really near to the 
border between the two areas, so one needs to be careful. 

Eventually, to find the second critical temperature, one 
computes the only solution of the equation in T c2 


T c2 -C 


N-M 


2 qh 2 


= 0 (C5) 


M T (qh 2 ,T c2 ) + 3g3 (e 7r B T <=2 ) 


if qh 2 < 77 , or 


T c2 -C 


M 


-2 qh 2 

53 ( 1 ) — g3(e k B T c2 ) 


= 0 (C6) 


if qh 2 > 77 . The bisection method allows to find this value 
with the requested accuracy, taking 0 for the lower bound 
and the temperature T c 1 for the upper bound. 
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